function[WCF2,WSQ2, Weffi] = welfare_fun_add(side, CFvec, cf, flexible)

global foldermats  cols  homedealer  kappa quantityrobus day loyaltyinW keepall largetrade tradesize fej zerobenefit loyalty 

cols = [3,9,4,1,2,5,8,6,7];  
CF=CFvec(cf);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% LOAD SQ AND CF WITH RETAILERS  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

retailer=1;
if zerobenefit==0
    temp=['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_' , num2str(loyalty)]]; 
else 
   temp=['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_', num2str(loyalty), '_zerobenefit']]; 
end 
load(fullfile([foldermats,temp,'.mat']))

PISQ_R    = PiSQ;
EUSQ_R   = uBSQ;
  
rhojSQ_R  = rhojSQ;
sjSQ_R     = sjSQ;
if homedealer==0, sHjSQ_R   = sjSQ; else sHjSQ_R    = sHjSQ;end 

if loyaltyinW==0 || homedealer==0
    WSQ2_R   = WSQ;
    WshrinkSQ2_R = WshrinkSQ;
    WSQxij_R  = WxijSQ;
    WSQxijH_R = WxijHSQ;
    WSQvdj_R  = WvdjSQ;
else 
    WSQ2_R   = WSQwH;
    WshrinkSQ2_R = WshrinkSQwH;
    WSQxij_R  = WxijSQwH;
    WSQxijH_R = WxijHSQwH;
    WSQvdj_R  = WvdjSQwH;
end 

if CFvec(cf)==6||CFvec(cf)==10||CFvec(cf)==14 %in the counterfact with r=0, use the welfare without r
    WSQ2_R       = WSQ;
    WshrinkSQ2_R = WshrinkSQ;
    WSQxij_R     = WxijSQ;
    WSQxijH_R    = WxijHSQ;
    WSQvdj_R     = WvdjSQ;   
end 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% LOAD SQ WITH INSTITUTIONALS 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
retailer=0;
if zerobenefit==0
    temp=['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_', num2str(loyalty)]];  
else 
   temp=['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_', num2str(loyalty), '_zerobenefit']];  
end 
load(fullfile([foldermats,temp,'.mat']))
    
    %compute best dealer -- Note: this is no longer used. Only works with certain specifications of the model
    if zerobenefit==0
        input_file=fullfile([foldermats,['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_', num2str(loyalty)]],'.mat']) ;
    else
        input_file=fullfile([foldermats,['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_', num2str(loyalty), '_zerobenefit']],'.mat']) ;
    end
    [idx_best, dealer_best, Weffi] = allocative_efficient(input_file);

  
q_data_I  =q_data; 
etaEjSQ_I = etaEjSQ;
etaDjSQ_I = etaDjSQ;
        

PISQ_I    = PiSQ ;
EUSQ_I    = uBSQ +uESQ;
rhojSQ_I  = rhojSQ;
sjSQ_I    = sjSQ;
s0jSQ_I   = s0jSQ;
if homedealer==0, sHjSQ_I    = sjSQ; else sHjSQ_I    = sHjSQ;end 

if loyaltyinW==0 || homedealer==0
    WSQ2_I       = WSQ;
    WshrinkSQ2_I = WshrinkSQ;
    WSQxij_I     = WxijSQ;
    WSQxijH_I    = WxijHSQ;
    WSQvdj_I     = WvdjSQ;
else
    WSQ2_I       = WSQwH;
    WshrinkSQ2_I = WshrinkSQwH;
    WSQxij_I     = WxijSQwH;
    WSQxijH_I    = WxijHSQwH;
    WSQvdj_I     = WvdjSQwH;
end 

if CFvec(cf)==6||CFvec(cf)==10||CFvec(cf)==14 %in the counterfact with r=0, use the welfare without r
    WSQ2_I       = WSQ;
    WshrinkSQ2_I = WshrinkSQ;
    WSQxij_I     = WxijSQ;
    WSQxijH_I    = WxijHSQ;
    WSQvdj_I     = WvdjSQ;   
end 

%compute welfare:
PISQ     = kappa*PISQ_R+(1-kappa)*PISQ_I;
EUSQ     = kappa*EUSQ_R+(1-kappa)*EUSQ_I;

%number of active dealers (assuming same number per retail and institution - this is an approxmiation!) 
Jt=kappa*sum(~isnan(EUSQ_R),2)+(1-kappa)*sum(~isnan(EUSQ_I),2);
avgJt = mean(Jt(Jt~=0)); % expected mass of investors


WSQ    = nansum(PISQ+EUSQ,2);
WSQ_I  = nansum(PISQ_I+EUSQ_I ,2);
WSQ_R  = nansum(PISQ_R+EUSQ_R,2);
    
WSQ2        = (1-kappa)*WSQ2_I   + kappa*WSQ2_R;
WshrinkSQ2  = (1-kappa)*WshrinkSQ2_I   + kappa*WshrinkSQ2_R;
WSQxij      = (1-kappa)*WSQxij_I + kappa*WSQxij_R;
WSQxijH     = (1-kappa)*WSQxijH_I + kappa*WSQxijH_R;
WSQvdj      = (1-kappa)*WSQvdj_I + kappa*WSQvdj_R;
     

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% COUNTERFACTUALS FOR INSTITUTIONALS  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   

retailer=0;
if zerobenefit==0
    temp=['/counterfac_',[num2str(side),'_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize),'_' num2str(CFvec(cf)),'_',num2str(flexible), '_', num2str(fej), '_', num2str(loyalty)]];
else 
    temp=['/counterfac_',[num2str(side),'_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize),'_' num2str(CFvec(cf)),'_',num2str(flexible), '_', num2str(fej), '_', num2str(loyalty), '_zerobenefit']];
end 
    load(fullfile([foldermats,temp,'.mat']))
    

PICF_I      = PiSQ ;
EUCF_I      = uBSQ +uESQ;
rhojcf21_I  = rhojSQ;
sjCF_I      = sjSQ;
s0jCF_I     = s0jSQ;
if homedealer==0, sHjcf21_I    = sjSQ; else sHjcf21_I    = sHjSQ;end 


if loyaltyinW==0 || homedealer==0
    WCF2_I       = WSQ;
    WshrinkCF2_I = WshrinkSQ;
    WCFxij_I     = WxijSQ;
    WCFxijH_I    = WxijHSQ;
    WCFvdj_I     = WvdjSQ;
else
    WCF2_I       = WSQwH;
    WshrinkCF2_I = WshrinkSQwH;
    WCFxij_I     = WxijSQwH;
    WCFxijH_I    = WxijHSQwH;
    WCFvdj_I     = WvdjSQwH;
end 


%Compute welfare:
PICF     = kappa*PISQ_R+(1-kappa)*PICF_I;
EUCF     = kappa*EUSQ_R+(1-kappa)*EUCF_I;

WCF     = nansum(PISQ+EUSQ,2);
WCF_I  = nansum(PICF_I+EUCF_I ,2);
WCF_R  = nansum(PISQ_R+EUSQ_R,2);
    
WCF2        = (1-kappa)*WCF2_I   + kappa*WSQ2_R;
WshrinkCF2  = (1-kappa)*WshrinkCF2_I   + kappa*WshrinkSQ2_R;
WCFxij      = (1-kappa)*WCFxij_I + kappa*WSQxij_R;
WCFxijH     = (1-kappa)*WCFxijH_I + kappa*WSQxijH_R;
WCFvdj      = (1-kappa)*WCFvdj_I + kappa*WSQvdj_R;
     
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% DISPLAY FINDINGS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   

%Welfare decomposition 
if homedealer==0
    WDcomp_I    =(WCFxij_I-WSQxij_I)./(WCFxij_I + WCFvdj_I - (WSQxij_I+WSQvdj_I))*100;
    WDcomp_R    =(WCFxij_R-WSQxij_R)./(WCFxij_R + WCFvdj_R - (WSQxij_R+WSQvdj_R))*100;
    WDcomp      =(WCFxij-WSQxij)./(WCFxij + WCFvdj - (WSQxij+WSQvdj))*100;
    WDcompv     =(WCFvdj-WSQvdj)./(WCFxij + WCFvdj - (WSQxij+WSQvdj))*100;
    WDcomp_both =[WDcomp_I,WDcomp_R, WDcomp,WDcompv];
else 
    DW      = WCF2    - WSQ2;       
    DW_xij  = WCFxij  - WSQxij;     %basequality
    DW_xijH = WCFxijH - WSQxijH;    %homedealer benefit 
    DW_vdj  = WCFvdj  - WSQvdj;     %dealer value

    %checking:
    DW_I      = WCF2_I    - WSQ2_I;       
    DW_xij_I  = WCFxij_I  - WSQxij_I;     %basequality
    DW_xijH_I = WCFxijH_I - WSQxijH_I;    %homedealer benefit 
    DW_vdj_I  = WCFvdj_I  - WSQvdj_I;     %dealer value

    DW_R      = WSQ2_R    - WSQ2_R;       
    DW_xij_R  = WSQxij_R  - WSQxij_R;     %basequality
    DW_xijH_R = WSQxijH_R - WSQxijH_R;    %homedealer benefit 
    DW_vdj_R  = WSQvdj_R  - WSQvdj_R;     %dealer value

    lines=find(~isnan(DW));
    
    mean([DW- (DW_xij+ DW_xijH+DW_vdj)])
    %--> should be zero
    
    % hist([WCF2_I - WSQ2_I- (WCFxij_I - WSQxij_I+ WCFxijH_I - WSQxijH_I+WCFvdj_I-WSQvdj_I)])
    %--> correct
    
end 

    %Change in quote
    disp('----------------------------------------------------------------')
    disp('QUOTE')
    disp('----------------------------------------------------------------')
    disp(['Increase in median quote:', num2str(nanmedian(nanmedian(q_data_I(:,cols)-q_data_I(:,cols))))]);
    
    %Change in expected surplus and profit
    disp('----------------------------------------------------------------')
    disp('PROFIT AND UTILITY')
    disp('----------------------------------------------------------------')
    disp(['Expected profit in SQ per trade for the dealer from retail investors is ', num2str(nanmedian(nanmean(PISQ_R,2)))])
    disp(['Expected profit in CF per trade for the dealer from retail investors is ', num2str(nanmedian(nanmean(PISQ_R,2)))])
    
    disp(['Expected surplus in SQ per trade for the retail investor is ', num2str(nanmedian(nanmean(EUSQ_R,2)))])
    disp(['Expected surplus in CF per trade for the retail investor is ', num2str(nanmedian(nanmean(EUSQ_R,2)))])

    disp(['Expected profit in SQ per trade for the dealer from institutional is ', num2str(nanmedian(nanmean(PISQ_I,2)))])
    disp(['Expected profit in CF per trade for the dealer from institutional is ', num2str(nanmedian(nanmean(PICF_I(PICF_I~=-Inf),2)))])
    
    disp(['Expected surplus in SQ per trade for the institutional investor is ', num2str(nanmedian(nanmean(EUSQ_I,2)))])
    disp(['Expected surplus in CF per trade for the institutional investor is ', num2str(nanmedian(nanmean(EUCF_I,2)))])


    %Change in welfare for average trade 
    disp('----------------------------------------------------------------')
    disp('WELFARE')
    disp('----------------------------------------------------------------')
    
    disp(['Total Expected welfare:  ', num2str(nanmean(WCF2)), ' bps'])
    disp(['Total Expected welfare:  ', num2str(nanmean(WshrinkCF2)), ' bps when shrinking dealer values '])

    disp(['Expected welfare changes by ', num2str(nanmean(WCF2 - WSQ2)), ' bps', 'out of ', num2str(nanmean( WSQ2)), ' bps'])
    disp(['Expected welfare changes by ', num2str(nanmean(WshrinkCF2 - WshrinkSQ2)), ' bps when shrinking dealer values out of ', num2str(nanmean(WshrinkSQ2)), ' bps'])


    diffW=nanmean((WCF2 - WSQ2)./WSQ2,2).*100; %average over dealers
    diffWshrink=nanmean((WshrinkCF2 - WshrinkSQ2)./WshrinkSQ2,2).*100; %average over dealers
    
    disp(['Expected welfare changes by ', num2str(nanmean(diffW(diffW~=-Inf))), '%.'])
    
    %Welfare decomposition with homedealer
    if homedealer==1
        disp([num2str(nanmean(DW_xijH(lines)./(DW_vdj(lines)+DW_xijH(lines)+DW_xij(lines))*100)), ' % of the welfare change comes from the homedealer benefit'])
        disp([num2str(nanmean(DW_vdj(lines)./(DW_vdj(lines)+DW_xijH(lines)+DW_xij(lines))*100)), ' % of the welfare change comes from matching to  dealers with higher values'])    
        disp([num2str(nanmean(DW_xij(lines)./(DW_vdj(lines)+DW_xijH(lines)+DW_xij(lines))*100)), ' % of the welfare change comes from matching to  dealers with better basequality values'])    
    end 
    
    disp(['Expected welfare changes by ', num2str(nanmean(diffWshrink(diffWshrink~=-Inf))), '% when shrinking dealer values.'])

    %Platform participation
    disp('----------------------------------------------------------------')
    disp('TRADING')
    disp('----------------------------------------------------------------')
    disp(['The probability that a retailer enters the platform is ', num2str(nanmedian(nanmean(100-rhojSQ_R*100,2))), ' in the SQ'])
    disp(['The probability that a retailer enters the platform is ', num2str(nanmedian(nanmean(100-rhojSQ_R*100,2))), ' in the CF'])
    disp(['The probability that an institutional enters the platform is ', num2str(nanmedian(nanmean(100-rhojSQ_I*100,2))), ' in the SQ'])
    disp(['The probability that an institutional enters the platform is ', num2str(nanmedian(nanmean(100-rhojcf21_I*100,2))), ' in the CF'])

    disp(['Probability to trade with the home dealer in SQ for I dealer is :  ', num2str(nanmean(nanmean(rhojSQ_I + sHjSQ_I))*100), '%'])
    disp(['Probability to trade with the home dealer in CF for I dealer is :  ', num2str(nanmean(nanmean(rhojcf21_I + sHjcf21_I))*100), '%'])
    disp(['Probability to trade with the home dealer in SQ for R dealer is : ', num2str(nanmean(nanmean(rhojSQ_R+sHjSQ_R))*100), '%'])
    disp(['Probability to trade with the home dealer in CF for R dealer is : ', num2str(nanmean(nanmean(rhojSQ_R+sHjSQ_R))*100), '%'])


    disp('----------------------------------------------------------------')
    disp('Welfare table')
    disp('----------------------------------------------------------------')
    if loyaltyinW==1
        output = (nanmean(WCF2)-62.0007)/62.0007*100;
        outputhrink = (nanmean(WshrinkCF2)-57.4838)/57.4838*100;
    else 
        output = (nanmean(WCF2)-36.6132)/36.6132*100;
        outputhrink = (nanmean(WshrinkCF2)-32.0963)/32.0963*100;
    end 
    disp(['Welfare change relative to decentralized market: ', num2str(output), '%'])
    disp(['Welfare change relative to decentralized market when shrinkng vDs: ', num2str(outputhrink), '%'])

